Análisis de Datos I
Tarea 11

Liberías

library(readr)
library(kableExtra)
library(dplyr)
library(caret)
library(traineR)

Pregunta 1

Complete las demostraciones de los Teoremas 2 y 4 de la presentación de la clase.

Teorema 2

Teorema 4

Pregunta 2

Diseñe un algoritmo en pseudocódigo para el Método del Análisis Discriminante Lineal según la teoría vista en clase.

Entrada:
- X: Matriz de datos de tamaño n por m

Salida:
- resultado: Vector que cuenta con:
-- G: matriz de centros de gravedad de cada clase
-- u: Los factores discriminates
-- lambda : valores propios

Paso Inicial: Separar los datos en dos bases
- X1 = datos con m_1 variables númericas
- X2 = código disyuntivo completo de las modalidades de la variable catégorica

Funcion ADL(X):

Paso 1: Calcular la matriz diagonal de cantidad de individuos en cada 
- clase D_G <- diag(n_g)

Paso 2: Calcular la matriz de centros de gravedades
- G = (D_G)^-1*(X2)^t*X1

Paso 3: Calcular la matriz de la cual se obtienen los factores discriminantes
y valores propios asociados
- FD_matriz <- (X1^t*X1)^-1*G^t*D_G*G

Paso 4: Obtener los vectores y valores propios de FD_matriz
- u <- eigen(FD_matriz)$vectors
- lambda <- eigen(FD_matriz)$values

Paso 5: Devolver matriz de centros de gravedad de cada clase,los factores discriminates
y valores propios
- devolver lista(G, u, lambda)

Fin Funció

Ejercicio 3

Programa en R el algoritmo disenando en el ejercicio anterior

Carga de datos

Inicialmente se procede cargando los datos que se utilizarán en el Análisis Discriminante Lineal.

Ejemplo_AD <- read.csv("Ejemplo_AD.csv", sep = ";", dec='.', header=T, 
                       row.names = 1, stringsAsFactors = T)

Tumores <- read.csv("tumores.csv", sep = ",", dec='.', header=T, 
                       row.names = 1, stringsAsFactors = T)

Partición de los datos

A continuación se procede con la partición de los datos en 2 sets, uno de entrenamiento con el 70% de los datos, y otro con un 30% para las pruebas.

# Se fija la semilla para reproducibilidad.
set.seed(2901)

# Se obtiene los indices para el primer set de train.
indices <- createDataPartition(Ejemplo_AD[,6], p = .70,list = F)

# Se obtiene la muestra para aprendizaje y prueba.
Ejemplo_AD_train <- Ejemplo_AD[indices,]
Ejemplo_AD_test <- Ejemplo_AD[-indices,]

Ahora se repite el proceso para el segundo set de datos.

# Se fija la semilla para reproducibilidad.
set.seed(2901)

# Se obtiene los indices para el primer set de train.
indices <- createDataPartition(Tumores[,17], p = .70,list = F)

# Se obtiene la muestra para aprendizaje y prueba.
Tumores_train <- Tumores[indices,]
Tumores_test <- Tumores[-indices,]

Creación del algoritmo con base en el pseudo codigo.

Ahora se procede con la creación de una función llamada ADL la cual recibe un dataframe (datos) con la información de entrenamiento.

ADL <- function(datos){
  # Se convierte la información recibida en un matriz para mejorar la forma en 
  # que se trabaja con los datos.
  X <- data.matrix(datos)
  
  # Se obtiene el indice la la variable categórica que está al final.
  indice <- ncol(datos)
  
  # Se divide la información.
  X1 <- X[,-indice]
  
  # Se obtiene las modalidades de la variable categórica.
  modalidades <- as.vector(unique(datos[,indice]))
  
  # Se obtiene el numero de modalidades que será útil.
  n_mod <- length(modalidades)
  
  # Se inicializa la matriz que contendrá los códigos disyuntivos.
  X2 <- matrix(, nrow = nrow(X), ncol = 0)
  
  # Se recorre cada modalidad y se crean los códigos disyuntivos respectivos.
  for(i in 1:n_mod){
    X2 <- cbind(X2, as.numeric(datos[indice] == modalidades[i]))
  }

  # Se crea la matriz diagonal D_g.
  D_g <- table(datos[,indice], datos[,indice])
  
  # Se obtiene la matriz inversa.
  Inv_D_g <- solve(D_g)
  
  # Se obtiene la matriz G con los centros de gravedad.
  G <- Inv_D_g %*% t(X2) %*% X1
  
  # Se crea la matriz de la cual se obtendrán los factores discriminantes.
  FD_matriz <- solve(t(X1) %*% X1) %*% t(G) %*% D_g %*% G
  
  
  # Se obtiene los valores y vectores propios de la matriz anterior.
  eigen <- eigen(FD_matriz)
  
  # Se crea una lista que contendrá toda la información del modelo para la 
  # clásificación de nuevos individuos.
  resultado <- list(Re(eigen[[1]][1:(n_mod-1)]), Re(eigen[[2]][,1:(n_mod-1)]),G)
  names(resultado) <- c("Valores propios", "Factores dicriminantes", 
                        "Matriz de CGs")
  
  return(resultado)
}

Pruebas con la base Ejemplo_AD

resultado <- ADL(Ejemplo_AD_train)
resultado
## $`Valores propios`
## [1] 0.9952515 0.7653630
## 
## $`Factores dicriminantes`
##             [,1]        [,2]
## [1,] -0.27231057  0.15635992
## [2,] -0.04373598 -0.07289682
## [3,] -0.47381610  0.31289403
## [4,]  0.82218860 -0.93270597
## [5,] -0.15309580 -0.04893383
## 
## $`Matriz de CGs`
##         RT1      RT2       RT3       RT4      RT5
## A 11.142857 6.185714 2.9428571 0.2857143 20.58571
## B  7.228571 5.485714 1.8285714 0.4285714 20.97143
## C  4.000000 6.242857 0.9428571 1.0285714 21.12857

Ahora se procede a obtener las coordenadas y a graficar.

z1 <- data.matrix(Ejemplo_AD_train[,-6]) %*% resultado[[2]][,1]
z2 <- data.matrix(Ejemplo_AD_train[,-6]) %*% resultado[[2]][,2]
coordenadas <- data.frame(z1, z2)
head(coordenadas)
##           z1         z2
## A1 -7.467266  0.3419385
## B1 -5.890088 -1.6350419
## A2 -7.103186  1.3406749
## B2 -5.937315  0.4688522
## C2 -3.462590 -0.9231316
## C3 -4.460751 -2.2514080
plot(coordenadas, col = as.numeric(Ejemplo_AD_train[,6]))

# Se genera el modelo.
modelo <- train.lda(VC~., data = Ejemplo_AD_train)
modelo
## Call:
## lda(VC ~ ., data = Ejemplo_AD_train)
## 
## Prior probabilities of groups:
##         A         B         C 
## 0.3333333 0.3333333 0.3333333 
## 
## Group means:
##         RT1      RT2       RT3       RT4      RT5
## A 11.142857 6.185714 2.9428571 0.2857143 20.58571
## B  7.228571 5.485714 1.8285714 0.4285714 20.97143
## C  4.000000 6.242857 0.9428571 1.0285714 21.12857
## 
## Coefficients of linear discriminants:
##             LD1         LD2
## RT1 -0.64978726  0.34348261
## RT2 -0.03092065  0.20736978
## RT3 -1.15680956 -0.71918814
## RT4  2.30118918  1.34704442
## RT5 -0.26127867 -0.02082875
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9955 0.0045
plot(modelo, col = as.numeric(Ejemplo_AD_train[,6]))

Pruebas con la base Tumores

resultado <- ADL(Tumores_train)
resultado
## $`Valores propios`
## [1] 13.87335
## 
## $`Factores dicriminantes`
##  [1] -1.594503e-04 -1.858376e-07  6.503070e-05  7.699501e-01 -2.116669e-05
##  [6]  7.380679e-08 -1.100928e-05 -3.955864e-01 -4.907293e-01  5.068602e-02
## [11]  2.165606e-03 -1.253470e-03  4.219568e-06  6.471504e-02  5.569606e-02
## [16] -2.740991e-03
## 
## $`Matriz de CGs`
##        media   varianza desviacion.estandar   entropia  asimetria    kurtosis
## 0 56.7331769 7066.90761          264.232242 13.5106678 248.426747 10614.19353
## 1  0.4586836   37.03272            1.272828  0.0597919   1.104296    40.24215
##     contraste     energia         asm homogeneidad disiminitud correlacion
## 0 1048.696670 13.85704768 13.39262412  14.03807333  9.68436072 12.48428200
## 1    5.518784  0.06343335  0.05863896   0.06591815  0.06731729  0.05764201
##         psnr        ssim         mse       dc
## 0 987.788123 13.56432460 0.579330101 6.003325
## 1   5.003442  0.06235764 0.005727337 0.000000

Ahora se procede a obtener las coordenadas y a graficar.

z1 <- data.matrix(Tumores_train[,-17]) %*% resultado[[2]]
coordenadas <- data.frame(z1)
head(coordenadas)
##                   z1
## Image1  -0.003422729
## Image2  -0.003264609
## Image3  -0.003233921
## Image9  -0.002942110
## Image10 -0.002846295
## Image12 -0.002786929
modelo <- train.lda(tipo~., data = Tumores_train)
modelo
## Call:
## lda(tipo ~ ., data = Tumores_train)
## 
## Prior probabilities of groups:
##          0          1 
## 0.06494961 0.93505039 
## 
## Group means:
##      media varianza desviacion.estandar  entropia asimetria kurtosis contraste
## 0 6.603462 533.1435            18.32433 0.8607971  15.89806 579.3482  79.45146
## 1 3.940748 490.8750            18.35386 0.9384656  17.25599 737.2733  72.84360
##     energia       asm homogeneidad disiminitud correlacion     psnr      ssim
## 0 0.9132214 0.8441988    0.9489940   0.9691368   0.8298462 72.03232 0.8977350
## 1 0.9625255 0.9302661    0.9750997   0.6726861   0.8671717 68.61283 0.9421926
##          mse        dc
## 0 0.08245390 0.0000000
## 1 0.04024089 0.4169974
## 
## Coefficients of linear discriminants:
##                               LD1
## media                4.186789e-01
## varianza             4.865830e-04
## desviacion.estandar -1.702792e-01
## entropia            -8.071862e+02
## asimetria            5.240808e-02
## kurtosis            -1.828464e-04
## contraste            2.756603e-02
## energia              3.979922e+02
## asm                  5.582098e+02
## homogeneidad        -1.288854e+02
## disiminitud         -5.439523e+00
## correlacion          3.209882e+00
## psnr                -1.008248e-02
## ssim                -1.617448e+02
## mse                 -1.391308e+02
## dc                   6.739918e+00

Conclusiones

Finalmente pese a que se fijo la semilla y se corrieron las mismas bases de datos, se puede notar que los resultados obtenidos variaron en gran medida, aunque en el gráfico del ejemplo con la base de datos Ejemplo_AD se nota que el procedimiento fue exitoso al separar a los individuos de test. Se puede teorizar que la diferencia en los resultados radica en el tratamiento de los datos que se da en la librería, pero según lo visto en la tería suministrada por el profesor, la metodología empleada en el algoritmo manual es también correcta.

Pregunta 4

Con las definiciones anteriores pruebe lo siguiente: Si V;VB;VW son las matrices de covarianza total, inter-clase intra-clase, respectivamente, entonces:

1. V = VB + VW

2. \(\sum_{s=1}^rq_sg_s = 0\). Es decir, \(rang(C_g) \leq r-1\)

3. \(rang(C_g) = rang(V_B)\)

Además, para la tabla de datos Ejemplo AD.csv calcule: gA, gB, gC, V, VB, VW y verifique que V = VB + VW

Primeramente, se carga la base de datos

Ejemplo_AD <-read.csv("Ejemplo_AD.csv",sep = ";",dec='.',header=T,row.names = 1,stringsAsFactors = T)

Las columnas de variables númericas deben estar centradas, por lo que, la media de cada una debe ser cero. Verifiquemos si lo están:

colMeans(Ejemplo_AD[,-6])
##       RT1       RT2       RT3       RT4       RT5 
##  7.213333  5.746667  1.950000  0.710000 21.740000

Se puede observar que las medias de cada columna númerica son diferente de cero, por lo que, debemos centrarlos de la siguiente manera:

for(i in 1: (ncol(Ejemplo_AD)-1)) {
  Ejemplo_AD[,i] <- Ejemplo_AD[,i]-mean(Ejemplo_AD[,i])
}

head(Ejemplo_AD)
##           RT1        RT2   RT3   RT4   RT5 VC
## A1  1.7866667 -1.1466667  0.05 -0.61  4.06  A
## B1 -3.2133333 -2.4466667 -1.55 -0.11 10.66  B
## C1 -5.8133333 -4.7466667 -0.85 -0.21  1.76  C
## A2  2.7866667  0.9533333  1.95 -0.51 -6.04  A
## B2  0.9866667  0.4533333  0.15 -0.51 -4.74  B
## C2 -4.7133333 -0.3466667 -1.15 -0.31 -5.44  C

Se obtienen los dos juegos de datos con las que se trabajan, la primera con variables númericas y la segunda con el código disyuntivo completo de las modalidades de la variable categórica.

X <- data.matrix(Ejemplo_AD)

# Datos con las variables númericas
X1<- X[,-6]
head(X1)
##           RT1        RT2   RT3   RT4   RT5
## A1  1.7866667 -1.1466667  0.05 -0.61  4.06
## B1 -3.2133333 -2.4466667 -1.55 -0.11 10.66
## C1 -5.8133333 -4.7466667 -0.85 -0.21  1.76
## A2  2.7866667  0.9533333  1.95 -0.51 -6.04
## B2  0.9866667  0.4533333  0.15 -0.51 -4.74
## C2 -4.7133333 -0.3466667 -1.15 -0.31 -5.44

El segundo juego de datos es:

# Datos con las variables categóricas
VC.A <- as.numeric(Ejemplo_AD$VC == "A")
VC.B <- as.numeric(Ejemplo_AD$VC == "B")
VC.C <- as.numeric(Ejemplo_AD$VC == "C")

X2 <- cbind(VC.A,VC.B,VC.C)
head(X2)
##      VC.A VC.B VC.C
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
## [4,]    1    0    0
## [5,]    0    1    0
## [6,]    0    0    1

Luego, se calcula matriz \(D_g\) diagonal del número de individuos en cada modalidad.

D_G <- table(Ejemplo_AD$VC,Ejemplo_AD$VC)
D_G
##    
##      A  B  C
##   A 10  0  0
##   B  0 10  0
##   C  0  0 10

Ahora, se obtienen los centros de gravedad de A, B y C.

\(g_A\)

# Matriz cuyas filas son los centros de gravedad
G <- solve(D_G) %*% t(X2) %*% X1

gA <- G[1,]
gA
##        RT1        RT2        RT3        RT4        RT5 
##  3.6866667  0.8433333  1.3400000 -0.3500000 -0.4000000

\(g_B\)

gB <- G[2,]
gB
##         RT1         RT2         RT3         RT4         RT5 
## -0.51333333  0.09333333 -0.21000000  0.25000000  0.48000000

\(g_C\)

gC <- G[3,]
gC
##        RT1        RT2        RT3        RT4        RT5 
## -3.1733333 -0.9366667 -1.1300000  0.1000000 -0.0800000

Seguido, se procede a calcular la matriz de covarianzas total \(V\).

# Matriz de pesos del conjuntos de individuos de la matriz X1
peso <-1/nrow(X1)
pesos_ind <- diag(peso, nrow(X1))

V <- t(X1)%*%pesos_ind %*%X1
kable_styling(kable(V))
RT1 RT2 RT3 RT4 RT5
RT1 10.8924889 1.1097111 3.328667 -0.6114667 -3.081200
RT2 1.1097111 10.8398222 2.405000 0.4715333 -4.521533
RT3 3.3286667 2.4050000 1.915833 -0.1175000 -2.416667
RT4 -0.6114667 0.4715333 -0.117500 0.5929000 1.017933
RT5 -3.0812000 -4.5215333 -2.416667 1.0179333 17.267067

La matriz de covarianzas inter-clase \(V_B\) es:

# Matriz diagonal de pesos de las 3 clases
pesos_clases <- D_G*peso

VB <- t(G)%*%pesos_clases%*%G
kable_styling(kable(VB))
RT1 RT2 RT3 RT4 RT5
RT1 7.9750222 2.0111778 2.8779333 -0.5786667 -0.4890667
RT2 2.0111778 0.5324222 0.7229667 -0.1218333 -0.0725333
RT3 2.8779333 0.7229667 1.0388667 -0.2115000 -0.1821333
RT4 -0.5786667 -0.1218333 -0.2115000 0.0650000 0.0840000
RT5 -0.4890667 -0.0725333 -0.1821333 0.0840000 0.1322667

La matriz de covarianzas intra-clase \(V_W\) es:

clases <- unique(Ejemplo_AD$VC)

# Se calcula la matriz VW
VW <- matrix(0, nrow = 5, ncol = 5)

for(i in (1:length(clases))) {
  clase <- Ejemplo_AD[Ejemplo_AD$VC == clases[i], -6]
  n <- nrow(clase)

  # Calcular la suma para la clase actual
  suma <- matrix(0, nrow = 5, ncol = 5)
  for(j in 1: n){
    xi_minus_g <- as.numeric(clase[j,] - G[i,])
    suma <- suma + xi_minus_g%*%t(xi_minus_g)
  }
  VW <- VW + suma
}
VW <- peso*VW
kable_styling(kable(VW))
2.9174667 -0.9014667 0.4507333 -0.0328000 -2.5921333
-0.9014667 10.3074000 1.6820333 0.5933667 -4.4490000
0.4507333 1.6820333 0.8769667 0.0940000 -2.2345333
-0.0328000 0.5933667 0.0940000 0.5279000 0.9339333
-2.5921333 -4.4490000 -2.2345333 0.9339333 17.1348000

Verificamos que se cumple \(V = V_B+V_W\).

kable_styling(kable(VB+VW))
RT1 RT2 RT3 RT4 RT5
RT1 10.8924889 1.1097111 3.328667 -0.6114667 -3.081200
RT2 1.1097111 10.8398222 2.405000 0.4715333 -4.521533
RT3 3.3286667 2.4050000 1.915833 -0.1175000 -2.416667
RT4 -0.6114667 0.4715333 -0.117500 0.5929000 1.017933
RT5 -3.0812000 -4.5215333 -2.416667 1.0179333 17.267067

Se puede notar que la suma de esas matrices es igual a \(V\).